function [data] = DGP(method,config,alternative)

N = config.N;
T = config.T;
n = config.n;
p = config.p;
spatial_model = config.spatial_model;
spatial_corr = config.spatial_corr;
serial_corr = config.serial_corr;
stationary = config.stationary;
discrete = config.discrete;
rho_endo = config.rho_endo;
pi_z = config.pi_z;

alpha = alternative;
beta = ones(p+1,1);

if N == 205
    pre_load = load('preload_data_small.mat');
else
    pre_load = load('preload_data_large.mat');
end
dis_mat = pre_load.dis_mat;


switch lower(spatial_model)
    case 'dampedcos'
        
        switch lower(method)
            case 'ols'
                theta = spatial_corr; % spatial correlation
                rho = serial_corr; % inter-temporal correlation
                rho_regressor = .5;

                % spatial and inter-temporal covariance 
                t_vec = kron(ones(N,1),[1;2]);
                Sigma = exp(-dis_mat/theta-squareform(pdist(t_vec))/rho);

                % covariance between regressors
                Sigma_control = (rho_regressor*ones(1+p)+(1-...
                    rho_regressor)*eye(1+p))/2; 

                regressor_matrix = chol(Sigma)*randn(n,1+p)*...
                    chol(Sigma_control)';
                D_star = regressor_matrix(:,1);
                
                if discrete == 1                    
                    c0 = 53.66/100;
                    q0 = quantile(D_star,c0);
                    c1 = c0+26.10/100;
                    q1 = quantile(D_star,c1);
                    c2 = c1+10.73/100;
                    q2 = quantile(D_star,c2);
                    c3 = c2+2.93/100;
                    q3 = quantile(D_star,c3);
                    c4 = c3+1.71/100;
                    q4 = quantile(D_star,c4);
                    c5 = c4+0.98/100;
                    q5 = quantile(D_star,c5);
                    c6 = c5+0.73/100;
                    q6 = quantile(D_star,c6);
                    c7 = c6+1.71/100;
                    q7 = quantile(D_star,c7);
                    c8 = c7+0.49/100;
                    q8 = quantile(D_star,c8);
                    c9 = c8+0.49/100;
                    q9 = quantile(D_star,c9);
                    f_D = @(x)(x>q0)+(x>q1)+(x>q2)+(x>q3)+(x>q4)+(x>q5)...
                        +(x>q6)+(x>q7)+(x>q8)+(x>q9);
                    D = f_D(D_star);
                    D = D/std(D)*std(D_star);
                else
                    D = D_star;
                end
                    
                if stationary == 1
                    X = [regressor_matrix(:,2:end),ones(n,1)];
                else
                    % location heterogeneity
                    coord = pre_load.coord;
                    temp = coord-repmat(mean(coord),length(coord),1);
                    hetero = 2*(sqrt(sum(temp.^2,2))-4);
                    X = [regressor_matrix(:,2:end).*...
                        repmat(hetero,1,p),ones(n,1)];
                end
                
                U = chol(Sigma)*randn(n,1);
                Y = D*alpha+X*beta+U;
             
                data.D = D;
                data.X = X;
                data.Y = Y;
                
            case 'iv'
                
                theta = spatial_corr; % spatial correlation
                rho = serial_corr; % inter-temporal correlation
                rho_regressor = .5;
%                 pi_z = .3; % relevance
                
                % spatial and inter-temporal covariance 
                t_vec = kron(ones(N,1),[1;2]);
                Sigma = exp(-dis_mat/theta-squareform(pdist(t_vec))/rho);
                
                % covariance between regressors
                Sigma_control = (rho_regressor*ones(1+p)+(1-...
                    rho_regressor)*eye(1+p))/2; 

                Z_mat = chol(Sigma)*randn(n,1+p)*...
                    chol(Sigma_control)';
     
                Z = Z_mat(:,1); % IV
                X = [Z_mat(:,2:end),ones(n,1)]; % controls
               
                Sigma_V = [eye(n) rho_endo*eye(n); rho_endo*eye(n) eye(n)];
                V = chol(Sigma_V)'*randn(2*n,1);
                V_Y = V(1:n);
                V_D = V((n+1):end);
                U_Y = chol(Sigma)'*V_Y;
                U_D = chol(Sigma)'*V_D;
                
                D_star = Z*pi_z+X*ones(p+1,1)+U_D;
                
                % tranform continuous to discrete
                % binary
%                 c = quantile(D_star,53.66/100);
%                 D = D_star>c;

                c0 = 53.66/100;
                q0 = quantile(D_star,c0);
                c1 = c0+26.10/100;
                q1 = quantile(D_star,c1);
                c2 = c1+10.73/100;
                q2 = quantile(D_star,c2);
                c3 = c2+2.93/100;
                q3 = quantile(D_star,c3);
                c4 = c3+1.71/100;
                q4 = quantile(D_star,c4);
                c5 = c4+0.98/100;
                q5 = quantile(D_star,c5);
                c6 = c5+0.73/100;
                q6 = quantile(D_star,c6);
                c7 = c6+1.71/100;
                q7 = quantile(D_star,c7);
                c8 = c7+0.49/100;
                q8 = quantile(D_star,c8);
                c9 = c8+0.49/100;
                q9 = quantile(D_star,c9);
                f_D = @(x)(x>q0)+(x>q1)+(x>q2)+(x>q3)+(x>q4)+(x>q5)...
                    +(x>q6)+(x>q7)+(x>q8)+(x>q9);
                D = f_D(D_star);
                D = D/std(D)*std(D_star);
%                 tabulate(D)

                Y = D*alpha+X*ones(p+1,1)+U_Y;
                
                data.D = D;
                data.X = X;
                data.Y = Y;
                data.Z = Z;
                
        end

    case 'sar' % SAR model
        
        dis_mat_unique = pre_load.dis_mat_unique;
        
        switch lower(method)
            case 'ols'                
                
                % SAR weighting matrix
%                 load distance_matrix
%                 d_bar = max(max(dis_mat_unique))/10;
                d_bar = .3;
                W = dis_mat_unique<d_bar;
                W = W-diag(diag(W));
                
                theta = spatial_corr;
                rho = exp(-serial_corr);
                rho_regressor = .5;
                
                % U
                V1 = randn(N,1);
                V2 = randn(N,1);
                U1 = (eye(N)-theta*W)\V1;
                U2 = rho*U1+sqrt(1-rho^2)*((eye(N)-theta*W)\V2);
                U = zeros(N*T,1);
                U(1:2:N*T) = U1;
                U(2:2:N*T) = U2;
                
                % covariance between regressors
                Sigma_control = (rho_regressor*ones(1+p)+(1-...
                    rho_regressor)*eye(1+p))/2; 
                
                % regressors
                V1 = randn(N,1+p)*chol(Sigma_control)';
                V2 = randn(N,1+p)*chol(Sigma_control)';
                X_mat = zeros(N*T,p+1);
                for i = 1 : p+1
                    X1 = (eye(N)-theta*W)\V1(:,i);
                    X2 = (eye(N)-theta*W)\V2(:,i);
                    X_mat(1:2:N*T,i) = X1;
                    X_mat(2:2:N*T,i) = X2;
                end
                
                D = X_mat(:,1);
                X = [X_mat(:,2:end),ones(n,1)];
                Y = D*alpha+X*beta+U;
                
                data.D = D;
                data.X = X;
                data.Y = Y;
                
            case 'iv'
        end
        
end






